# Analysis of Study 1 in 
# Voters Use Campaign Finance Transparency and Compliance Information
# Abby K. Wood
# Political Behavior, 2022

# This file contains the code to replicate all findings and figures
# for Study 1 in the manuscript (the vignette experiment with candidates 
# Johnson and Conley), as well as robustness analysis in 
# Appendices 2 (CCES), 4, 5, and 10A.

# it uses the following datasets, downloadable on the Political Behavior dataverse
# figure1_df.RData
# full_replication_data_s1.RData
# no_extremes_replication_data_s1.RData

# set working directory
# setwd()


# install packages
#install.packages("ri2")
#install.packages("car")
#install.packages("grDevices")
#install.packages("ggplot2")

# load packages
library(car)
library(grDevices)
library(ri2)
library(ggplot2)
library(reshape)


# FIGURE 1 (RELATIVE IMPORTANCE OF VALENCE ISSUES) ####
# load CCES15 item CFT_361_1 - CFT_361_6, which I've cleaned into the 
# relative importance rannkings, all of which are coded 1 = very important,
# 2 = important,  3 = moderately important, 4 = of little importance, and
# 5 = unimportant.  6 = NA (did not answer)

# imp.cmoney = importance of amount campaign has raised (CFT361_1)
# imp.gmoney = importance of amount outside groups have spent 
##              in support of candidate (CFT361_2)
# imp.disc = importance of the amount candidate discloses re financing (CFT361_3)
# imp.pers = importance of persuasiveness in public (CFT361_4)
# imp.grasp  = importance of grasp of the issues (CFT361_5)
# imp.prof  = imprtance of professional background (x$CFT361_6)

load("figure1_df.RData")
df$r <- as.numeric(1:2000)
df1 <- melt(df, id.vars = 'r', variable_name = 'item', na.rm=F)
df1$item <- as.factor(df1$item)
df1 <- subset(df1, !(is.na(value)))

# collapse values in melted dataset to "Important or Very Important", 
# "Moderately Important", "Unimportant or Of Little Importance"
df1$value <- ifelse(df1$value < 3, "Unimpt / Little Imptce", ifelse(
        df1$value == 3, "Moderately Impt", ifelse(
                df1$value > 3, "Impt / Very Impt", NA)))
levels(df1$value) <- c("Impt / Very Impt",
                       "Moderately Impt", 
                       "Unimpt / Little Imptce")
df1$response <- df1$value

# black and white color scheme
bwPalette <- c( "#525252", "#969696", "#cccccc")

g <- ggplot(df1, aes(x=factor(item), fill = response))
g + geom_bar(position ="fill") +
        theme_bw() +
        scale_fill_manual(values=bwPalette) +
        coord_cartesian(ylim = c(0, 1)) +
        labs(title = "", 
             x = "", y = "Percent") +
        scale_x_discrete(labels = c("camp. money", "group money", "disclosure", 
                                    "persuasiveness", "prof background", 
                                    "grasp of issues")) +
        theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = 'black', 
                                         size = 13))
#dev.off()

# remove those data 
rm(df, df1, g, bwPalette)


## STUDY 1 REPLICATION BEGINS HERE ####

load("full_replication_data_s1.RData")
load("no_extremes_replication_data_s1.RData")


# STUDY 1 CODEBOOK ####

# jfav = johnson favorability
# cfav = conley favorability
# netfav = johnson - conley favorability
# jvote = likelihood of voting for johnson
# cvote = likelihood of voting for conley
# netvote = jvote - cvote
# jtrust = johnson trustworthiness
# ctrust = conley trustworthiness
# nettrust = jtrust - ctrust

# tr takes four values:  
## 1 = conley +/- 10, control
## 2 = conley +/- 20, control
## 3 = conley +/- 10, treatment
## 4 = conley +/- 20, treatment

# race
## 1 = white
## 2 = black
## 3 = hispanic
## 4 = asian
## 5 = native american
## 6 = mixed
## 7 = other
## 8 = middle eastern
## 98 = skipped
## 99 = not asked

# CC15_350 = party
## 1 = independent
## 2 = democratic
## 3 = republican

# USC450 = self-identified ideology 
## (1-10 in response to how conservative/liberal one is,
## where "conservative" or "liberal" is asked based on 
## one's party ID or ideology, and only true independents
## with party ID of independent AND 4 on 7 point ideology scale
## are randomized into either the "how conservative" or "how liberal"
## version of the question.)

# CC15_340a = 7 point ideology, where
## 1 = extremely liberal,
## 2 = liberal
## 3 = somewhat liberal
## 4 = moderate
## 5 = somewhat conservative
## 6 = conservative
## 7 = extremely conservative, 

# CC15_310a congressional control question 1

# CC15_310b congressional control question 2

# creating and re-naming a few variables
x$ideo7pt <- x$CC15_340a
x$ideo7pt[x$ideo7pt > 7] <- NA
# this is an indicator of whether one is extremely liberal or 
# extremely conservative based on the 7 point ideology score
x$ext <- as.numeric(x$ideo7 == 1 | x$ideo7 == 7)

# this is the 1-10 ranking of conservatism or liberalism, with true
# independents randomly assigned to either conservatism or liberalism
x$extreme <- x$USC450

# creating two party dataset
x$republican <- as.numeric(x$CC15_350 == 3)
x$democrat <- as.numeric(x$CC15_350 == 2)
x$independent <- as.numeric(x$CC15_350 == 1)
dems <- subset(x, democrat == 1)
reps <- subset(x, republican == 1)
inds <- subset(x, independent == 1)
twop <- rbind(dems, reps)
twop$ideo2p <- NA
twop$ideo2p[twop$democrat == 1] <- -1*twop$extreme[twop$democrat == 1]
twop$ideo2p[twop$republican == 1] <- twop$extreme[twop$republican == 1]

twop <- twop[twop$tr == 2 | twop$tr == 4,]
twop$treat <- as.numeric(twop$tr == 4)


# Study 1 estimates in text surrounding Figure 3 ####
# as I analyze, I create a data.frame called "plotdata" that has columns for
# (1) effect type (e.g., "Johnson Favorability"), 
# (2) an estimate (diff in means), 
# (3) ci.low (from t.test), 
# (4) ci.high (from t.test)


# "Respondents told that Johnson has dark money support disapproved 
# (45.9 and 45.6 on favorability and vote choice)." 

jf <- t.test(x$jfav[x$tr == 4], x$jfav[x$tr == 2])
jv <- t.test(x$jvote[x$tr == 4], x$jvote[x$tr == 2])
jf
jv

# "The resulting favorability gap is -9.45 (95% CI = [-13.03, -5.90]); 
# vote choice gap is -9.19 [-13.03, -5.35]."
jf$estimate[1] - jf$estimate[2]
jf$conf.int

jv$estimate[1] - jv$estimate[2]
jv$conf.int

# create plotdata matrix
plotdata <- matrix(nrow = 6, ncol = 8)
colnames(plotdata) <- c("effect", "estimate", "ci.low", "ci.high", "pval", 
                        "mean_tr", "mean_ct", "arrow_plot_name")

# create plotdata vector for johnson favorability for Figure 3
plotdata[1,1] <- "Johnson - Favorability (T-C)"
plotdata[1,2] <- jf$estimate[1] - jf$estimate[2]
plotdata[1,3] <- jf$conf.int[1]
plotdata[1,4] <- jf$conf.int[2]
plotdata[1,5] <- round(jf$p.value, digits = 2)
plotdata[1,6] <- jf$estimate[1]
plotdata[1,7] <- jf$estimate[2]
plotdata[1,8] <- "Johnson - Favorability"

# create johnson vote vector for Figure 3
plotdata[2,1] <- "Johnson - Probable Vote (T-C)"
plotdata[2,2] <- jv$estimate[1] - jv$estimate[2]
plotdata[2,3] <- jv$conf.int[1]
plotdata[2,4] <- jv$conf.int[2]
plotdata[2,5] <- round(jv$p.value, digits = 2)
plotdata[2,6] <- jv$estimate[1]
plotdata[2,7] <- jv$estimate[2]
plotdata[2,8] <- "Johnson - Probable Vote"


# "Respondents in the treatment group rated Conley 59.9 for favorability, a 
# difference of 4.6 [0.88, 8.32] from the control group."
cf <- t.test(x$cfav[x$tr == 4], x$cfav[x$tr == 2])
cf$estimate[1] - cf$estimate[2]
cf$conf.int

# "Vote choice results are similar, at 4.15 [0.32, 7.99]." 
cv <- t.test(x$cvote[x$tr == 4], x$cvote[x$tr == 2])
cv$estimate[1] - cv$estimate[2]
cv$conf.int

# conley favorability plot vector
plotdata[3,1] <- "Conley - Favorability (T-C)"
plotdata[3,2] <- cf$estimate[1] - cf$estimate[2]
plotdata[3,3] <- cf$conf.int[1]
plotdata[3,4] <- cf$conf.int[2]
plotdata[3,5] <- round(cf$p.value, digits = 2)
plotdata[3,6] <- cf$estimate[1]
plotdata[3,7] <- cf$estimate[2]
plotdata[3,8] <- "Conley - Favorability"

# conley vote plot vector
plotdata[4,1] <- "Conley - Probable Vote (T-C)"
plotdata[4,2] <- cv$estimate[1] - cv$estimate[2]
plotdata[4,3] <- cv$conf.int[1]
plotdata[4,4] <- cv$conf.int[2]
plotdata[4,5] <- round(cv$p.value, digits = 2)
plotdata[4,6] <- cv$estimate[1]
plotdata[4,7] <- cv$estimate[2]
plotdata[4,8] <- "Conley - Probable Vote"

# "Net favorability is estimated to be -13.99 [-18.47, -9.47],"
nf <- t.test(x$netfav[x$tr == 4], x$netfav[x$tr == 2])
nf$estimate[1] - nf$estimate[2]
nf$conf.int

# "...and net vote likelihood -13.3 [-18.24, -8.43]." 
nv <- t.test(x$netvote[x$tr == 4], x$netvote[x$tr == 2])
nv$estimate[1] - nv$estimate[2]
nv$conv.int

# net favorability plotdata vector
plotdata[5,1] <- "Net Favorability (J-C)"
plotdata[5,2] <- nf$estimate[1] - nf$estimate[2]
plotdata[5,3] <- nf$conf.int[1]
plotdata[5,4] <- nf$conf.int[2]
plotdata[5,5] <- round(nf$p.value, digits = 2)
plotdata[5,6] <- nf$estimate[1]
plotdata[5,7] <- nf$estimate[2]
plotdata[5,8] <- "Net Favorability (J-C)"

# net vote plotdata vector
plotdata[6,1] <- "Net Vote (J-C)"
plotdata[6,2] <- nv$estimate[1] - nv$estimate[2]
plotdata[6,3] <- nv$conf.int[1]
plotdata[6,4] <- nv$conf.int[2]
plotdata[6,5] <- round(nv$p.value, digits = 2)
plotdata[6,6] <- nv$estimate[1]
plotdata[6,7] <- nv$estimate[2]
plotdata[6,8] <- "Net Vote (J-C)"

# put all plot data together as a dataframe
plotdata <- as.data.frame(plotdata)
# hard code p values that are zero, all of which are at least e-06
plotdata.tmp <- rep(NA, length(plotdata$pval))
plotdata.tmp <- ifelse(plotdata$pval == "0", "0.0000", as.character(plotdata$pval))
plotdata$pval <- plotdata.tmp
plotdata$mean_tr <- as.numeric(as.character(plotdata$mean_tr))
plotdata$mean_ct <- as.numeric(as.character(plotdata$mean_ct))


# Footnote: The control group was equally favorable between Johnson and Conley (55.3) 
# and was slightly more likely to vote for Conley (58.4) over Johnson 
# (54.8, difference not statistically significant).
fn1 <- t.test(x$jfav[x$tr == 2], x$cfav[x$tr == 2])
fn1$estimate[1]
fn1$estimate[2]

fn2 <- t.test(x$jvote[x$tr == 2], x$cvote[x$tr == 2])
fn2$estimate[2]
fn2$estimate[1]
fn2$conf.int


# Figure 3 ####

plot(x = NULL, y = NULL, axes = F, 
     # set values of different parameters
     xlim = c(-48, 20),
     ylim = c(1, 6.5),
     xlab = "",
     ylab = "",
     main = "")

        # set margins and letter size
        par(cex=0.8, mai = c(0.5, 0.35, 0.5, 0.35))

        # add the 0 vertical line and tick marks
        abline(v = 0)
        #segments(x0 = 0, y0 = 0, x1 = 0, y1 = 6.5)
        axis(side=1,at=c(-20, -10, 0, 10, 20),tick=TRUE, las=2, cex.axis=0.7)

        # Fill the figure with the information which is inside the 'plotdata' dataframe
        # Add the estimates as points
        points(as.character(plotdata$estimate[6]), 1, type = 'p')
        points(as.character(plotdata$estimate[5]), 2, type = 'p')
        points(as.character(plotdata$estimate[4]), 3, type = 'p')
        points(as.character(plotdata$estimate[3]), 4, type = 'p')
        points(as.character(plotdata$estimate[2]), 5, type = 'p')
        points(as.character(plotdata$estimate[1]), 6, type = 'p')


        # Add the confidence intervals as lines
        segments(x0=as.numeric(as.character(plotdata$ci.low[6])), y0=1, 
                 x1=as.numeric(as.character(plotdata$ci.high[6])), y1=1)
        segments(x0=as.numeric(as.character(plotdata$ci.low[5])), y0=2, 
                 x1=as.numeric(as.character(plotdata$ci.high[5])), y1=2)
        segments(x0=as.numeric(as.character(plotdata$ci.low[4])), y0=3, 
                 x1=as.numeric(as.character(plotdata$ci.high[4])), y1=3)
        segments(x0=as.numeric(as.character(plotdata$ci.low[3])), y0=4, 
                 x1=as.numeric(as.character(plotdata$ci.high[3])), y1=4)
        segments(x0=as.numeric(as.character(plotdata$ci.low[2])), y0=5, 
                 x1=as.numeric(as.character(plotdata$ci.high[2])), y1=5)
        segments(x0=as.numeric(as.character(plotdata$ci.low[1])), y0=6, 
                 x1=as.numeric(as.character(plotdata$ci.high[1])), y1=6)


        # Add each variable name
        text(-45, 6, plotdata[1,1], adj = 0, cex=1)
        text(-45, 5, plotdata[2,1], adj = 0, cex=1)
        text(-45, 4, plotdata[3,1], adj = 0, cex=1)
        text(-45, 3, plotdata[4,1], adj = 0, cex=1)
        text(-45, 2, plotdata[5,1], adj = 0, cex=1)
        text(-45, 1, plotdata[6,1], adj = 0, cex=1)

        # Add pvalues
        text(18, 6, plotdata$pval[1], cex=1)
        text(18, 5, plotdata$pval[2], cex=1)
        text(18, 4, plotdata$pval[3], cex=1)
        text(18, 3, plotdata$pval[4], cex=1)
        text(18, 2, plotdata$pval[5], cex=1)
        text(18, 1, plotdata$pval[6], cex=1)

        # Add horizontal lines to make it pretty
        abline(h = 2.5, lty = 3)
        abline(h = 4.5, lty = 3)

        # Add the subtitles  
        text(18, 6.5, "P Value", font = 2)
        text(-42, 6.5, "Outcome", font = 2)
        mtext("Percentage Point Change", side = 1, line=2, 
              cex = 0.7, at = c(0,0))






# STUDY 1 EXPLORATORY ANALYSIS AND FIGURE 4

# first, self reported ideology on 7 pt scale (x$ideo7pt)--
x$treatment <- as.numeric(x$tr == 4 | x$tr == 3)
x$treatment <- ifelse(x$treatment == 1, "treatment", "control")
x$id <- 1:length(x$treatment)

# subset just to the 20pp diffs
x20pp <- subset(x, tr == 2 | tr == 4)

# aggregate to create dataframes for each panel.
jfav <- aggregate(x20pp$jfav, by = list(x20pp$treatment, x20pp$ideo7pt), 
                  FUN = "mean", na.rm=T)
names(jfav) <- c("treatment", "ideo7pt", "jfav")
cfav <- aggregate(x20pp$cfav, by = list(x20pp$treatment, x20pp$ideo7pt), 
                  FUN = "mean", na.rm=T)
names(cfav) <- c("treatment", "ideo7pt", "cfav")

jcfav <- cbind(jfav, cfav$cfav)
jcfavct <- subset(jcfav, treatment == "control")
jcfavtr <- subset(jcfav, treatment == "treatment")

jvote <- aggregate(x20pp$jvote, by = list(x20pp$treatment, x20pp$ideo7pt), 
                   FUN = "mean", na.rm=T)
names(jvote) <- c("treatment", "ideo7pt", "jvote")
cvote <- aggregate(x20pp$cvote, by = list(x20pp$treatment, x20pp$ideo7pt), 
                   FUN = "mean", na.rm=T)
names(cvote) <- c("treatment", "ideo7pt", "cvote")

jcvote <- cbind(jvote, cvote$cvote)
jcvotect <- subset(jcvote, treatment == "control")
jcvotetr <- subset(jcvote, treatment == "treatment")


ctdata <- cbind(jcfavct, jcvotect$jvote, jcvotect$cvote)
names(ctdata) <- c("treatment", "ideo7pt", "jfav", "cfav", "jvote", "cvote")
trdata <- cbind(jcfavtr, jcvotetr$jvote, jcvotetr$cvote)
names(trdata) <- c("treatment", "ideo7pt", "jfav", "cfav", "jvote", "cvote")




# Figure 4 ####
# Figure 4 panels are created separately. 
# I saved as 4 separate .pdfs, FWIW, but I took the pdf() commands out so as not 
# to force anyone using this code to save files on their drives.

        # Johnson vote choice for Fig 4
        plot(x = NULL, y = NULL, axes = F, 
             # set values of different parameters
             xlim = c(1, 7),
             ylim = c(40, 80),
             ylab = "Vote Choice - Johnson",
             xlab = "Ideology",
             main = "")
        
        # add the 0 vertical line and tick marks
        abline(h = 0)
        #segments(x0 = 0, y0 = 0, x1 = 0, y1 = 6.5)
        axis(side=1,at=c(1:7),tick=TRUE, las=2, cex.axis=0.7)
        axis(side=2,at=c(40, 50, 60, 70, 80),tick=TRUE, las=2, cex.axis=0.7)
        # Fill the figure with the information which is inside the 'results' dataframe
        
        # First, add the estimates as arrows FOR JOHNSON
        arrows(x0 = 7, y0 = ctdata$jvote[7], 
               x1 = 7, y1 = trdata$jvote[7], length = 0.1)
        arrows(x0 = 6, y0 = ctdata$jvote[6], 
               x1 = 6, y1 = trdata$jvote[6], length = 0.1)
        arrows(x0 = 5, y0 = ctdata$jvote[5], 
               x1 = 5, y1 = trdata$jvote[5], length = 0.1)
        arrows(x0 = 4, y0 = ctdata$jvote[4], 
               x1 = 4, y1 = trdata$jvote[4], length = 0.1)
        arrows(x0 = 3, y0 = ctdata$jvote[3], 
               x1 = 3, y1 = trdata$jvote[3], length = 0.1)
        arrows(x0 = 2, y0 = ctdata$jvote[2], 
               x1 = 2, y1 = trdata$jvote[2], length = 0.1)
        arrows(x0 = 1, y0 = ctdata$jvote[1], 
               x1 = 1, y1 = trdata$jvote[1], length = 0.1)

        
        ## Conley vote chocie for Figure 4
        plot(x = NULL, y = NULL, axes = F, 
             # set values of different parameters
             xlim = c(1, 7),
             ylim = c(40, 80),
             ylab = "Vote Choice - Conley",
             xlab = "Ideology",
             main = "")
        
        # add the 0 vertical line and tick marks
        abline(h = 0)
        #segments(x0 = 0, y0 = 0, x1 = 0, y1 = 6.5)
        axis(side=1,at=c(1:7),tick=TRUE, las=2, cex.axis=0.7)
        axis(side=2,at=c(40, 50, 60, 70, 80),tick=TRUE, las=2, cex.axis=0.7)
        # Fill the figure with the information which is inside the 'results' dataframe
        
        # First, add the estimates as points FOR JOHNSON
        arrows(x0 = 7, y0 = ctdata$cvote[7], 
               x1 = 7, y1 = trdata$cvote[7], length = 0.1)
        arrows(x0 = 6, y0 = ctdata$cvote[6], 
               x1 = 6, y1 = trdata$cvote[6], length = 0.1)
        arrows(x0 = 5, y0 = ctdata$cvote[5], 
               x1 = 5, y1 = trdata$cvote[5], length = 0.1)
        arrows(x0 = 4, y0 = ctdata$cvote[4], 
               x1 = 4, y1 = trdata$cvote[4], length = 0.1)
        arrows(x0 = 3, y0 = ctdata$cvote[3], 
               x1 = 3, y1 = trdata$cvote[3], length = 0.1)
        arrows(x0 = 2, y0 = ctdata$cvote[2], 
               x1 = 2, y1 = trdata$cvote[2], length = 0.1)
        arrows(x0 = 1, y0 = ctdata$cvote[1], 
               x1 = 1, y1 = trdata$cvote[1], length = 0.1)
        
        
        ####
        # Johnson favorability for Fig 4
          plot(x = NULL, y = NULL, axes = F, 
             # set values of different parameters
             xlim = c(1, 7),
             ylim = c(40, 80),
             ylab = "Favorability - Johnson",
             xlab = "Ideology",
             main = "")
        
        # add the 0 vertical line and tick marks
        abline(h = 0)
        #segments(x0 = 0, y0 = 0, x1 = 0, y1 = 6.5)
        axis(side=1,at=c(1:7),tick=TRUE, las=2, cex.axis=0.7)
        axis(side=2,at=c(40, 50, 60, 70, 80),tick=TRUE, las=2, cex.axis=0.7)
        # Fill the figure with the information which is inside the 'results' dataframe
        
        # First, add the estimates as arrows FOR JOHNSON
        arrows(x0 = 7, y0 = ctdata$jfav[7], 
               x1 = 7, y1 = trdata$jfav[7], length = 0.1)
        arrows(x0 = 6, y0 = ctdata$jfav[6], 
               x1 = 6, y1 = trdata$jfav[6], length = 0.1)
        arrows(x0 = 5, y0 = ctdata$jfav[5], 
               x1 = 5, y1 = trdata$jfav[5], length = 0.1)
        arrows(x0 = 4, y0 = ctdata$jfav[4], 
               x1 = 4, y1 = trdata$jfav[4], length = 0.1)
        arrows(x0 = 3, y0 = ctdata$jfav[3], 
               x1 = 3, y1 = trdata$jfav[3], length = 0.1)
        arrows(x0 = 2, y0 = ctdata$jfav[2], 
               x1 = 2, y1 = trdata$jfav[2], length = 0.1)
        arrows(x0 = 1, y0 = ctdata$jfav[1], 
               x1 = 1, y1 = trdata$jfav[1], length = 0.1)
        
   
        # Conley favorability for Fig 4
        plot(x = NULL, y = NULL, axes = F, 
             # set values of different parameters
             xlim = c(1, 7),
             ylim = c(40, 80),
             ylab = "Favorability - Conley",
             xlab = "Ideology",
             main = "")
        
        # add the 0 vertical line and tick marks
        abline(h = 0)
        #segments(x0 = 0, y0 = 0, x1 = 0, y1 = 6.5)
        axis(side=1,at=c(1:7),tick=TRUE, las=2, cex.axis=0.7)
        axis(side=2,at=c(40, 50, 60, 70, 80),tick=TRUE, las=2, cex.axis=0.7)
        # Fill the figure with the information which is inside the 'results' dataframe
        
        # First, add the estimates as arrows FOR JOHNSON
        arrows(x0 = 7, y0 = ctdata$cfav[7], 
               x1 = 7, y1 = trdata$cfav[7], length = 0.1)
        arrows(x0 = 6, y0 = ctdata$cfav[6], 
               x1 = 6, y1 = trdata$cfav[6], length = 0.1)
        arrows(x0 = 5, y0 = ctdata$cfav[5], 
               x1 = 5, y1 = trdata$cfav[5], length = 0.1)
        arrows(x0 = 4, y0 = ctdata$cfav[4], 
               x1 = 4, y1 = trdata$cfav[4], length = 0.1)
        arrows(x0 = 3, y0 = ctdata$cfav[3], 
               x1 = 3, y1 = trdata$cfav[3], length = 0.1)
        arrows(x0 = 2, y0 = ctdata$cfav[2], 
               x1 = 2, y1 = trdata$cfav[2], length = 0.1)
        arrows(x0 = 1, y0 = ctdata$cfav[1], 
               x1 = 1, y1 = trdata$cfav[1], length = 0.1)
        
  
# Estimates in text surrounding Figure 4 refer the reader to Appendix 5c, 
        # so please refer to the replication code for Appendix 5c 
        # to replicate those.
        

# APPENDICES RELATED TO STUDY 1 ####

# APPENDIX 2, PANEL DEMOGRAPHICS WITH STUDY 1 SAMPLE. ####
        
# APPENDIX 4A - REPLICATION OF STUDY 1, OMITTING MORE MODERATE CONLEY ####
 
        # use dataframe y, loaded above.  It's the smaller dataset omitting 
        # cases in which Respondent was more extreme than Conley
        
        # Johnson favorability
        jf <- t.test(y$jfav[y$tr == 4], y$jfav[y$tr == 2])
        
        plotdata <- matrix(nrow = 6, ncol = 8)
        colnames(plotdata) <- c("effect", "estimate", "ci.low", "ci.high", "pval", 
                                "mean_tr", "mean_ct", "arrow_plot_name")
        plotdata[1,1] <- "Johnson - Favorability (T-C)"
        plotdata[1,2] <- jf$estimate[1] - jf$estimate[2]
        plotdata[1,3] <- jf$conf.int[1]
        plotdata[1,4] <- jf$conf.int[2]
        plotdata[1,5] <- round(jf$p.value, digits = 2)
        plotdata[1,6] <- jf$estimate[1]
        plotdata[1,7] <- jf$estimate[2]
        plotdata[1,8] <- "Johnson - Favorability"
        
        # johnson vote
        jv <- t.test(y$jvote[y$tr == 4], y$jvote[y$tr == 2])
        
        plotdata[2,1] <- "Johnson - Probable Vote (T-C)"
        plotdata[2,2] <- jv$estimate[1] - jv$estimate[2]
        plotdata[2,3] <- jv$conf.int[1]
        plotdata[2,4] <- jv$conf.int[2]
        plotdata[2,5] <- round(jv$p.value, digits = 2)
        plotdata[2,6] <- jv$estimate[1]
        plotdata[2,7] <- jv$estimate[2]
        plotdata[2,8] <- "Johnson - Probable Vote"
        
        # Conley favorability
        cf <- t.test(y$cfav[y$tr == 4], y$cfav[y$tr == 2])
        
        plotdata[3,1] <- "Conley - Favorability (T-C)"
        plotdata[3,2] <- cf$estimate[1] - cf$estimate[2]
        plotdata[3,3] <- cf$conf.int[1]
        plotdata[3,4] <- cf$conf.int[2]
        plotdata[3,5] <- round(cf$p.value, digits = 2)
        plotdata[3,6] <- cf$estimate[1]
        plotdata[3,7] <- cf$estimate[2]
        plotdata[3,8] <- "Conley - Favorability"
        
        # Conley vote
       cv <- t.test(y$cvote[y$tr == 4], y$cvote[y$tr == 2])
        
        plotdata[4,1] <- "Conley - Probable Vote (T-C)"
        plotdata[4,2] <- cv$estimate[1] - cv$estimate[2]
        plotdata[4,3] <- cv$conf.int[1]
        plotdata[4,4] <- cv$conf.int[2]
        plotdata[4,5] <- round(cv$p.value, digits = 2)
        plotdata[4,6] <- cv$estimate[1]
        plotdata[4,7] <- cv$estimate[2]
        plotdata[4,8] <- "Conley - Probable Vote"
        
        
        # net favorability
        nf <- t.test(y$netfav[y$tr == 3], y$netfav[y$tr == 1])
        
        plotdata[5,1] <- "Net Favorability (J-C)"
        plotdata[5,2] <- nf$estimate[1] - nf$estimate[2]
        plotdata[5,3] <- nf$conf.int[1]
        plotdata[5,4] <- nf$conf.int[2]
        plotdata[5,5] <- round(nf$p.value, digits = 2)
        plotdata[5,6] <- nf$estimate[1]
        plotdata[5,7] <- nf$estimate[2]
        plotdata[5,8] <- "Net Favorability (J-C)"
        
        
        # net vote
        nv <- t.test(y$netvote[y$tr == 3], y$netvote[y$tr == 1])
        
        plotdata[6,1] <- "Net Vote (J-C)"
        plotdata[6,2] <- nv$estimate[1] - nv$estimate[2]
        plotdata[6,3] <- nv$conf.int[1]
        plotdata[6,4] <- nv$conf.int[2]
        plotdata[6,5] <- round(nv$p.value, digits = 2)
        plotdata[6,6] <- nv$estimate[1]
        plotdata[6,7] <- nv$estimate[2]
        plotdata[6,8] <- "Net Vote (J-C)"
        
        
        plotdata <- as.data.frame(plotdata)
        # hard code p values that are zero, all of which are at least e-06
        plotdata.tmp <- rep(NA, length(plotdata$pval))
        plotdata.tmp <- ifelse(plotdata$pval == "0", "0.0000", as.character(plotdata$pval))
        plotdata$pval <- plotdata.tmp
        plotdata$mean_tr <- as.numeric(as.character(plotdata$mean_tr))
        plotdata$mean_ct <- as.numeric(as.character(plotdata$mean_ct))
        
        
       # plot A4A
        plot(x = NULL, y = NULL, axes = F, 
             # set values of different parameters
             xlim = c(-48, 20),
             ylim = c(1, 6.5),
             xlab = "",
             ylab = "",
             main = "")
        # set margins and letter size
        par(cex=0.8, mai = c(0.5, 0.35, 0.5, 0.35))
        
        # add the 0 vertical line and tick marks
        abline(v = 0)
        #segments(x0 = 0, y0 = 0, x1 = 0, y1 = 6.5)
        axis(side=1,at=c(-20, -10, 0, 10, 20),tick=TRUE, las=2, cex.axis=0.7)
        
        # add x axis label
        mtext("Percentage Points", side = 1, line=2, cex = 0.7, at = c(0,0))
        # Fill the figure with the information which is inside the 'results' dataframe
        
        # add the estimates as points
        points(as.character(plotdata$estimate[6]), 1, type = 'p')
        points(as.character(plotdata$estimate[5]), 2, type = 'p')
        points(as.character(plotdata$estimate[4]), 3, type = 'p')
        points(as.character(plotdata$estimate[3]), 4, type = 'p')
        points(as.character(plotdata$estimate[2]), 5, type = 'p')
        points(as.character(plotdata$estimate[1]), 6, type = 'p')
        
        
        # add the confidence intervals as lines
        segments(x0=as.numeric(as.character(plotdata$ci.low[6])), y0=1, 
                 x1=as.numeric(as.character(plotdata$ci.high[6])), y1=1)
        segments(x0=as.numeric(as.character(plotdata$ci.low[5])), y0=2, 
                 x1=as.numeric(as.character(plotdata$ci.high[5])), y1=2)
        segments(x0=as.numeric(as.character(plotdata$ci.low[4])), y0=3, 
                 x1=as.numeric(as.character(plotdata$ci.high[4])), y1=3)
        segments(x0=as.numeric(as.character(plotdata$ci.low[3])), y0=4, 
                 x1=as.numeric(as.character(plotdata$ci.high[3])), y1=4)
        segments(x0=as.numeric(as.character(plotdata$ci.low[2])), y0=5, 
                 x1=as.numeric(as.character(plotdata$ci.high[2])), y1=5)
        segments(x0=as.numeric(as.character(plotdata$ci.low[1])), y0=6, 
                 x1=as.numeric(as.character(plotdata$ci.high[1])), y1=6)
        
        
        # add each variable name
        text(-45, 6, plotdata[1,1], adj = 0, cex=1)
        text(-45, 5, plotdata[2,1], adj = 0, cex=1)
        text(-45, 4, plotdata[3,1], adj = 0, cex=1)
        text(-45, 3, plotdata[4,1], adj = 0, cex=1)
        text(-45, 2, plotdata[5,1], adj = 0, cex=1)
        text(-45, 1, plotdata[6,1], adj = 0, cex=1)
        
        # add pvalues
        text(18, 6, plotdata$pval[1], cex=1)
        text(18, 5, plotdata$pval[2], cex=1)
        text(18, 4, plotdata$pval[3], cex=1)
        text(18, 3, plotdata$pval[4], cex=1)
        text(18, 2, plotdata$pval[5], cex=1)
        text(18, 1, plotdata$pval[6], cex=1)
        
        # add horizontal lines to make it pretty
        abline(h = 2.5, lty = 3)
        abline(h = 4.5, lty = 3)
        
        # add the subtitles  
        text(18, 6.5, "P Value", font = 2)
        text(-42, 6.5, "Outcome", font = 2)
        
      
        
# APPENDIX 4B - REPLICATION OF STUDY 1, 10PP DIFFERENCE B/W CANDIDATES ####

        # as before, creating a data.frame called <plotdata> that has columns for
        # (1) effect type (e.g., "Johnson Favorability"), 
        # (2) an estimate (diff in means), 
        # (3) ci.low (from t.test), 
        # (4) ci.high (from t.test)
        
        # johnson favorability results, 10pp separation
        jf <- t.test(x$jfav[x$tr == 3], x$jfav[x$tr == 1])
        
        # creating dataframe and including these results
        plotdata <- matrix(nrow = 6, ncol = 8)
        colnames(plotdata) <- c("effect", "estimate", "ci.low", "ci.high", "pval", 
                                "mean_tr", "mean_ct", "arrow_plot_name")
        plotdata[1,1] <- "Johnson - Favorability (T-C)"
        plotdata[1,2] <- jf$estimate[1] - jf$estimate[2]
        plotdata[1,3] <- jf$conf.int[1]
        plotdata[1,4] <- jf$conf.int[2]
        plotdata[1,5] <- round(jf$p.value, digits = 2)
        plotdata[1,6] <- jf$estimate[1]
        plotdata[1,7] <- jf$estimate[2]
        plotdata[1,8] <- "Johnson - Favorability"
        
        # johnson vote results, 10pp separation
        jv <- t.test(x$jvote[x$tr == 3], x$jvote[x$tr == 1])
        
        # including johnson vote results, 10pp separation, in plot data
        plotdata[2,1] <- "Johnson - Probable Vote (T-C)"
        plotdata[2,2] <- jv$estimate[1] - jv$estimate[2]
        plotdata[2,3] <- jv$conf.int[1]
        plotdata[2,4] <- jv$conf.int[2]
        plotdata[2,5] <- round(jv$p.value, digits = 2)
        plotdata[2,6] <- jv$estimate[1]
        plotdata[2,7] <- jv$estimate[2]
        plotdata[2,8] <- "Johnson - Probable Vote"
        
        # conley favorability, 10pp separation
        cf <- t.test(x$cfav[x$tr == 3], x$cfav[x$tr == 1])
        
        # including conley favorability vector, 10pp separation, in plot data
        plotdata[3,1] <- "Conley - Favorability (T-C)"
        plotdata[3,2] <- cf$estimate[1] - cf$estimate[2]
        plotdata[3,3] <- cf$conf.int[1]
        plotdata[3,4] <- cf$conf.int[2]
        plotdata[3,5] <- round(cf$p.value, digits = 2)
        plotdata[3,6] <- cf$estimate[1]
        plotdata[3,7] <- cf$estimate[2]
        plotdata[3,8] <- "Conley - Favorability"
        
        # conley vote, 10pp separation
        cv <- t.test(x$cvote[x$tr == 3], x$cvote[x$tr == 1])
        
        # including conley vote, 10pp separation, in plot
        plotdata[4,1] <- "Conley - Probable Vote (T-C)"
        plotdata[4,2] <- cv$estimate[1] - cv$estimate[2]
        plotdata[4,3] <- cv$conf.int[1]
        plotdata[4,4] <- cv$conf.int[2]
        plotdata[4,5] <- round(cv$p.value, digits = 2)
        plotdata[4,6] <- cv$estimate[1]
        plotdata[4,7] <- cv$estimate[2]
        plotdata[4,8] <- "Conley - Probable Vote"
        
        
        # net favorability, 10pp separation
        nf <- t.test(x$netfav[x$tr == 3], x$netfav[x$tr == 1])
        
        # including net favorability, 10pp separation, in plot
        plotdata[5,1] <- "Net Favorability (J-C)"
        plotdata[5,2] <- nf$estimate[1] - nf$estimate[2]
        plotdata[5,3] <- nf$conf.int[1]
        plotdata[5,4] <- nf$conf.int[2]
        plotdata[5,5] <- round(nf$p.value, digits = 2)
        plotdata[5,6] <- nf$estimate[1]
        plotdata[5,7] <- nf$estimate[2]
        plotdata[5,8] <- "Net Favorability (J-C)"
        
        # net vote, 10pp separation
        nv <- t.test(x$netvote[x$tr == 3], x$netvote[x$tr == 1])
        
        # including net vote, 10pp separation in plot
        plotdata[6,1] <- "Net Vote (J-C)"
        plotdata[6,2] <- nv$estimate[1] - nv$estimate[2]
        plotdata[6,3] <- nv$conf.int[1]
        plotdata[6,4] <- nv$conf.int[2]
        plotdata[6,5] <- round(nv$p.value, digits = 2)
        plotdata[6,6] <- nv$estimate[1]
        plotdata[6,7] <- nv$estimate[2]
        plotdata[6,8] <- "Net Vote (J-C)"
        
        # compile plot data
        plotdata <- as.data.frame(plotdata)
        # hard code p values that are zero, all of which are at least e-06
        plotdata.tmp <- rep(NA, length(plotdata$pval))
        plotdata.tmp <- ifelse(plotdata$pval == "0", "0.0000", as.character(plotdata$pval))
        plotdata$pval <- plotdata.tmp
        plotdata$mean_tr <- as.numeric(as.character(plotdata$mean_tr))
        plotdata$mean_ct <- as.numeric(as.character(plotdata$mean_ct))
        
        
        
    # Now make the actual plot in Appendix 4B ----
        
        plot(x = NULL, y = NULL, axes = F, 
             # set values of different parameters
             xlim = c(-48, 20),
             ylim = c(1, 6.5),
             xlab = "",
             ylab = "",
             main = "")
        # set margins and letter size
        par(cex=0.8, mai = c(0.5, 0.35, 0.5, 0.35))
        
        # add the 0 vertical line and tick marks
        abline(v = 0)
        #segments(x0 = 0, y0 = 0, x1 = 0, y1 = 6.5)
        axis(side=1,at=c(-20, -10, 0, 10, 20),tick=TRUE, las=2, cex.axis=0.7)
        # Fill the figure with the information which is inside the 'results' dataframe
        
        # First, add the estimates as points
        points(as.character(plotdata$estimate[6]), 1, type = 'p')
        points(as.character(plotdata$estimate[5]), 2, type = 'p')
        points(as.character(plotdata$estimate[4]), 3, type = 'p')
        points(as.character(plotdata$estimate[3]), 4, type = 'p')
        points(as.character(plotdata$estimate[2]), 5, type = 'p')
        points(as.character(plotdata$estimate[1]), 6, type = 'p')
        
        
        # Then, add the confidence intervals as lines
        segments(x0=as.numeric(as.character(plotdata$ci.low[6])), y0=1, 
                 x1=as.numeric(as.character(plotdata$ci.high[6])), y1=1)
        segments(x0=as.numeric(as.character(plotdata$ci.low[5])), y0=2, 
                 x1=as.numeric(as.character(plotdata$ci.high[5])), y1=2)
        segments(x0=as.numeric(as.character(plotdata$ci.low[4])), y0=3, 
                 x1=as.numeric(as.character(plotdata$ci.high[4])), y1=3)
        segments(x0=as.numeric(as.character(plotdata$ci.low[3])), y0=4, 
                 x1=as.numeric(as.character(plotdata$ci.high[3])), y1=4)
        segments(x0=as.numeric(as.character(plotdata$ci.low[2])), y0=5, 
                 x1=as.numeric(as.character(plotdata$ci.high[2])), y1=5)
        segments(x0=as.numeric(as.character(plotdata$ci.low[1])), y0=6, 
                 x1=as.numeric(as.character(plotdata$ci.high[1])), y1=6)
        
        
        # Third, add each variable name
        text(-45, 6, plotdata[1,1], adj = 0, cex=1)
        text(-45, 5, plotdata[2,1], adj = 0, cex=1)
        text(-45, 4, plotdata[3,1], adj = 0, cex=1)
        text(-45, 3, plotdata[4,1], adj = 0, cex=1)
        text(-45, 2, plotdata[5,1], adj = 0, cex=1)
        text(-45, 1, plotdata[6,1], adj = 0, cex=1)
        
        # fourth, add pvalues
        text(18, 6, plotdata$pval[1], cex=1)
        text(18, 5, plotdata$pval[2], cex=1)
        text(18, 4, plotdata$pval[3], cex=1)
        text(18, 3, plotdata$pval[4], cex=1)
        text(18, 2, plotdata$pval[5], cex=1)
        text(18, 1, plotdata$pval[6], cex=1)
        
        # add horizontal lines to make it pretty
        abline(h = 2.5, lty = 3)
        abline(h = 4.5, lty = 3)
        
        # finally, add the subtitles  
        text(18, 6.5, "P Value", font = 2)
        text(-42, 6.5, "Outcome", font = 2)
        
        mtext("Percentage Points", side = 1, line=2, cex = 0.7, at = c(0,0))
        
       
# APPENDIX 5A - PARTISAN DIFFERENCES IN REWARD/PUNISHMENT ####
        # H4a - heterogeneous treatment effects by party
        
        # Table A5A, Col 1
        
        # johnson favorability, republicans
        jfavreps <- t.test(twop$jfav[twop$republican == 1 & twop$tr == 4], 
               twop$jfav[twop$republican == 1 & twop$tr == 2])
        jfavreps$estimate[1] - jfavreps$estimate[2]
        jfavreps$p.value
        
        # johnson favorability, democrats
        jfavdems <- t.test(twop$jfav[twop$democrat == 1 & twop$tr == 4], 
                           twop$jfav[twop$democrat == 1 & twop$tr == 2])
        jfavdems$estimate[1] - jfavdems$estimate[2]
        jfavdems$p.value
        
        # johnson favorability, f-stat between republicans and democrats
        # testing this using code from #5 on 
        # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        
        set.seed(555)
        
        # first on jfav
        twopjfav <- subset(twop, !(is.na(jfav)))
        
        # note that jfav has 319 complete cells, hence 319 here.
        declaration <- declare_ra(319)
        sims <- 10000
        
        twopjfav$Z <- twopjfav$treat
        
        ate_obs <- with(twopjfav, mean(jfav[republican == 1]) - mean(jfav[republican == 0]))
        
        ri_out_jfav <-
                conduct_ri(
                        model_1 = jfav ~ Z + republican,
                        model_2 = jfav ~ Z + republican + Z * republican,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = twopjfav, 
                        sims = sims
                )
        # From main text:  "Republicans were more punitive toward Johnson 
        # for favorability (F-statistic 5.26, p = 0.02)"
        summary(ri_out_jfav)
    
        # create a vector for replication purposes
        TA5Ac1 <- c(
          (jfavreps$estimate[1] - jfavreps$estimate[2]), 
          jfavreps$p.value, 
          (jfavdems$estimate[1] - jfavdems$estimate[2]), 
          jfavdems$p.value, 
          summary(ri_out_jfav)[[1, 2]], 
          summary(ri_out_jfav)[[1, 3]])

        # Table A5A, Col 2
        
        # johnson vote, republicans
        jvotereps <- t.test(twop$jvote[twop$republican == 1 & twop$tr == 4], 
                           twop$jvote[twop$republican == 1 & twop$tr == 2])
        jvotereps$estimate[1] - jvotereps$estimate[2]
        jvotereps$p.value
        
        # johnson vote, democrats
        jvotedems <- t.test(twop$jvote[twop$democrat == 1 & twop$tr == 4], 
                           twop$jvote[twop$democrat == 1 & twop$tr == 2])
        jvotedems$estimate[1] - jvotedems$estimate[2]
        jvotedems$p.value
        
        # johnson vote, f-stat between republicans and democrats
        # testing this using code from #5 on 
        # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        
        set.seed(555)
        
        # first on jvote
        twopjvote <- subset(twop, !(is.na(jvote)))
        
        # note that jvote has 328 complete cells, hence 328 here.
        declaration <- declare_ra(328)
        sims <- 10000
        
        twopjvote$Z <- twopjvote$treat
        
        ate_obs <- with(twopjvote, mean(jvote[republican == 1]) - mean(jvote[republican == 0]))
        
        ri_out_jvote <-
                conduct_ri(
                        model_1 = jvote ~ Z + republican,
                        model_2 = jvote ~ Z + republican + Z * republican,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = twopjvote, 
                        sims = sims
                )
        summary(ri_out_jvote)

        # create a vector for replication purposes
        TA5Ac2 <- c(
          (jvotereps$estimate[1] - jvotereps$estimate[2]), 
          jvotereps$p.value, 
          (jvotedems$estimate[1] - jvotedems$estimate[2]), 
          jvotedems$p.value, 
          summary(ri_out_jvote)[[1, 2]], 
          summary(ri_out_jvote)[[1, 3]])
        
        # Table A5A, Col 3
        # conley favorability, republicans
        cfavreps <- t.test(twop$cfav[twop$republican == 1 & twop$tr == 4], 
                           twop$cfav[twop$republican == 1 & twop$tr == 2])
        cfavreps$estimate[1] - cfavreps$estimate[2]
        cfavreps$p.value
        
        # conley favorability, democrats
        cfavdems <- t.test(twop$cfav[twop$democrat == 1 & twop$tr == 4], 
                           twop$cfav[twop$democrat == 1 & twop$tr == 2])
        cfavdems$estimate[1] - cfavdems$estimate[2]
        cfavdems$p.value
        
        # conley favorability, f-stat between republicans and democrats
        # testing this using code from #5 on 
        # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        
        set.seed(555)
        
        # first on cfav
        twopcfav <- subset(twop, !(is.na(cfav)))
        
        # note that cfav has 321 complete cells, hence 321 here.
        declaration <- declare_ra(321)
        sims <- 10000
        
        twopcfav$Z <- twopcfav$treat
        
        ate_obs <- with(twopcfav, mean(cfav[republican == 1]) - mean(cfav[republican == 0]))
        
        ri_out_cfav <-
                conduct_ri(
                        model_1 = cfav ~ Z + republican,
                        model_2 = cfav ~ Z + republican + Z * republican,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = twopcfav, 
                        sims = sims
                )
        summary(ri_out_cfav)
        
        # create a vector for replication purposes
        TA5Ac3 <- c(
          (cfavreps$estimate[1] - cfavreps$estimate[2]), 
          cfavreps$p.value, 
          (cfavdems$estimate[1] - cfavdems$estimate[2]), 
          cfavdems$p.value, 
          summary(ri_out_cfav)[[1, 2]], 
          summary(ri_out_cfav)[[1, 3]])
        
        
        # Table A5A, Col 4
        
        # conley vote, republicans
        cvotereps <- t.test(twop$cvote[twop$republican == 1 & twop$tr == 4], 
                            twop$cvote[twop$republican == 1 & twop$tr == 2])
        cvotereps$estimate[1] - cvotereps$estimate[2]
        cvotereps$p.value
        
        # conley vote, democrats
        cvotedems <- t.test(twop$cvote[twop$democrat == 1 & twop$tr == 4], 
                            twop$cvote[twop$democrat == 1 & twop$tr == 2])
        cvotedems$estimate[1] - cvotedems$estimate[2]
        cvotedems$p.value
        
        # conley vote, f-stat between republicans and democrats
        # testing this using code from #5 on 
        # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        
        set.seed(555)
        
        # first on cvote
        twopcvote <- subset(twop, !(is.na(cvote)))
        
        # note that cvote has 328 complete cells, hence 328 here.
        declaration <- declare_ra(328)
        sims <- 10000
        
        twopcvote$Z <- twopcvote$treat
        
        ate_obs <- with(twopcvote, mean(cvote[republican == 1]) - mean(cvote[republican == 0]))
        
        ri_out_cvote <-
                conduct_ri(
                        model_1 = cvote ~ Z + republican,
                        model_2 = cvote ~ Z + republican + Z * republican,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = twopcvote, 
                        sims = sims
                )
        # From main text:  "...and Democratic respondents may be more likely 
        # to vote for Conley (F-statistic = 2.18, p = 0.14)"
        summary(ri_out_cvote)

        # create a vector for replication purposes
        TA5Ac4 <- c(
          (cvotereps$estimate[1] - cvotereps$estimate[2]), 
          cvotereps$p.value, 
          (cvotedems$estimate[1] - cvotedems$estimate[2]), 
          cvotedems$p.value, 
          summary(ri_out_cvote)[[1, 2]], 
          summary(ri_out_cvote)[[1, 3]])
        
    # Here's the table:
        TA5A <- cbind(TA5Ac1, TA5Ac2, TA5Ac3, TA5Ac4)
        colnames(TA5A) <- c("Johnson fav", "Johnson vote", "Conley fav", "Conley vote")
        TA5A <- round(TA5A, digits = 2)
        
# APPENDIX 5B - RANDOMIZATION INFERENCE ON EXTREMISTS VS. EVERYONE ELSE ####
       
# table A5B, col 1
        
         # johnson favorability, extremists
        jfavext <- t.test(x20pp$jfav[x20pp$ext == 1 & x20pp$tr == 4], 
                           x20pp$jfav[x20pp$ext == 1 & x20pp$tr == 2])
        jfavext$estimate[1] - jfavext$estimate[2]
        jfavext$p.value
        
        # johnson favorability, non-extremists
        jfavnoext <- t.test(x20pp$jfav[x20pp$ext == 0 & x20pp$tr == 4], 
                          x20pp$jfav[x20pp$ext == 0 & x20pp$tr == 2])
        jfavnoext$estimate[1] - jfavnoext$estimate[2]
        jfavnoext$p.value
        
        x20pp$Z <- as.numeric(x20pp$tr == 4)
        # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        #install.packages("ri2")
        #library(ri2)
        set.seed(1000)
        
        # first on jfav
        xjfav <- subset(x20pp, !(is.na(jfav)) & !(is.na(ext)))
        
        # note that jfav has 458 complete cells
        declaration <- declare_ra(458)
        sims <- 10000
        
        ate_obs <- with(xjfav, mean(jfav[ext == 1]) - mean(jfav[ext == 0]))
        
        ri_out_jfav_ext <-
                conduct_ri(
                        model_1 = jfav ~ Z + ext + Z * ext,
                        model_2 = jfav ~ Z + ext,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = xjfav, 
                        sims = sims
                )
        
        summary(ri_out_jfav_ext)  # this p value makes no sense, given the plot. 
        plot(ri_out_jfav_ext) # this is the plot
        # We know that the p value on an F stat of 7.47 with DF of 455 and 456 
        # is < 0.0001, so that's what I report in the table.  Can't explain the
        # problem in the conduct_ri command here.
        
        TA5Bc1 <- c(
          (jfavext$estimate[1] - jfavext$estimate[2]), 
          jfavext$p.value, 
          (jfavnoext$estimate[1] - jfavnoext$estimate[2]), 
          jfavnoext$p.value, 
          summary(ri_out_jfav_ext)[[1, 2]], 
          0.00)
        
        
# table A5B, col 2 
        
        # johnson vote, extremists
        jvoteext <- t.test(x20pp$jvote[x20pp$ext == 1 & x20pp$tr == 4], 
                          x20pp$jvote[x20pp$ext == 1 & x20pp$tr == 2])
        jvoteext$estimate[1] - jvoteext$estimate[2]
        jvoteext$p.value
        
        # johnson vote, non-extremists
        jvotenoext <- t.test(x20pp$jvote[x20pp$ext == 0 & x20pp$tr == 4], 
                            x20pp$jvote[x20pp$ext == 0 & x20pp$tr == 2])
        jvotenoext$estimate[1] - jvotenoext$estimate[2]
        jvotenoext$p.value
        
        x20pp$Z <- as.numeric(x20pp$tr == 4)
        # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        #install.packages("ri2")
        #library(ri2)
        set.seed(1000)
        
        # first on jvote
        xjvote <- subset(x20pp, !(is.na(jvote)) & !(is.na(ext)))
        
        # note that jvote has 475 complete cells
        declaration <- declare_ra(475)
        sims <- 10000
        
        ate_obs <- with(xjvote, mean(jvote[ext == 1]) - mean(jvote[ext == 0]))
        
        ri_out_jvote_ext <-
                conduct_ri(
                        model_1 = jvote ~ Z + ext + Z * ext,
                        model_2 = jvote ~ Z + ext,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = xjvote, 
                        sims = sims
                )
        
        summary(ri_out_jvote_ext)  
        plot(ri_out_jvote_ext) 
        
        TA5Bc2 <- c(
          (jvoteext$estimate[1] - jvoteext$estimate[2]), 
          jvoteext$p.value, 
          (jvotenoext$estimate[1] - jvotenoext$estimate[2]), 
          jvotenoext$p.value, 
          summary(ri_out_jvote_ext)[[1, 2]], 
          summary(ri_out_jvote_ext)[[1, 3]])
        
        
# table A5B, col 3
        
        # conley favorability, extremists
        cfavext <- t.test(x20pp$cfav[x20pp$ext == 1 & x20pp$tr == 4], 
                          x20pp$cfav[x20pp$ext == 1 & x20pp$tr == 2])
        cfavext$estimate[1] - cfavext$estimate[2]
        cfavext$p.value
        
        # conley favorability, non-extremists
        cfavnoext <- t.test(x20pp$cfav[x20pp$ext == 0 & x20pp$tr == 4], 
                            x20pp$cfav[x20pp$ext == 0 & x20pp$tr == 2])
        cfavnoext$estimate[1] - cfavnoext$estimate[2]
        cfavnoext$p.value
        
        x20pp$Z <- as.numeric(x20pp$tr == 4)
        # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        #install.packages("ri2")
        #library(ri2)
        set.seed(1000)
        
        # first on cfav
        xcfav <- subset(x20pp, !(is.na(cfav)) & !(is.na(ext)))
        
        # note that cfav has 461 complete cells
        declaration <- declare_ra(461)
        sims <- 10000
        
        ate_obs <- with(xcfav, mean(cfav[ext == 1]) - mean(cfav[ext == 0]))
        
        ri_out_cfav_ext <-
                conduct_ri(
                        model_1 = cfav ~ Z + ext + Z * ext,
                        model_2 = cfav ~ Z + ext,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = xcfav, 
                        sims = sims
                )
        
        summary(ri_out_cfav_ext)  
        
        
        TA5Bc3 <- c(
          (cfavext$estimate[1] - cfavext$estimate[2]), 
          cfavext$p.value, 
          (cfavnoext$estimate[1] - cfavnoext$estimate[2]), 
          cfavnoext$p.value, 
          summary(ri_out_cfav_ext)[[1, 2]], 
          summary(ri_out_cfav_ext)[[1, 3]])
        
        
# table A5B, col 4
        
        # conley vote, extremists
        cvoteext <- t.test(x20pp$cvote[x20pp$ext == 1 & x20pp$tr == 4], 
                           x20pp$cvote[x20pp$ext == 1 & x20pp$tr == 2])
        cvoteext$estimate[1] - cvoteext$estimate[2]
        cvoteext$p.value
        
        # conley vote, non-extremists
        cvotenoext <- t.test(x20pp$cvote[x20pp$ext == 0 & x20pp$tr == 4], 
                             x20pp$cvote[x20pp$ext == 0 & x20pp$tr == 2])
        cvotenoext$estimate[1] - cvotenoext$estimate[2]
        cvotenoext$p.value
        
      # https://egap.org/resource/10-randomization-inference-procedures-with-ri2/
        #install.packages("ri2")
        #library(ri2)
        set.seed(1000)
        
        # first on cvote
        xcvote <- subset(x20pp, !(is.na(cvote)) & !(is.na(ext)))
        
        # note that cvote has 475 complete cells
        declaration <- declare_ra(475)
        sims <- 10000
        
        ate_obs <- with(xcvote, mean(cvote[ext == 1]) - mean(cvote[ext == 0]))
        
        ri_out_cvote_ext <-
                conduct_ri(
                        model_1 = cvote ~ Z + ext + Z * ext,
                        model_2 = cvote ~ Z + ext,
                        declaration = declaration,
                        sharp_hypothesis = ate_obs,
                        data = xcvote, 
                        sims = sims
                )
        
        summary(ri_out_cvote_ext)  
        
        TA5Bc4 <- c(
          (cvoteext$estimate[1] - cvoteext$estimate[2]), 
          cvoteext$p.value, 
          (cvotenoext$estimate[1] - cvotenoext$estimate[2]), 
          cvotenoext$p.value, 
          summary(ri_out_cvote_ext)[[1, 2]], 
          summary(ri_out_cvote_ext)[[1, 3]])
        
        
        # Here's the table:
        TA5B <- cbind(TA5Bc1, TA5Bc2, TA5Bc3, TA5Bc4)
        colnames(TA5B) <- c("Johnson fav", "Johnson vote", "Conley fav", "Conley vote")
        TA5B <- round(TA5B, digits = 2)
        
        
# APPENDIX 5C - IDEOLOGICAL DIFFERENCES IN REWARD/PUNISHMENT ####
# (Table to support Fig 4)

        jv1 <- t.test(x$jvote[x$tr == 4 & x$ideo7pt == 1], x$jvote[x$tr == 2 & x$ideo7pt == 1])
        jv2 <- t.test(x$jvote[x$tr == 4 & x$ideo7pt == 2], x$jvote[x$tr == 2 & x$ideo7pt == 2])
        jv3 <- t.test(x$jvote[x$tr == 4 & x$ideo7pt == 3], x$jvote[x$tr == 2 & x$ideo7pt == 3])
        jv4 <- t.test(x$jvote[x$tr == 4 & x$ideo7pt == 4], x$jvote[x$tr == 2 & x$ideo7pt == 4])
        jv5 <- t.test(x$jvote[x$tr == 4 & x$ideo7pt == 5], x$jvote[x$tr == 2 & x$ideo7pt == 5])
        jv6 <- t.test(x$jvote[x$tr == 4 & x$ideo7pt == 6], x$jvote[x$tr == 2 & x$ideo7pt == 6])
        jv7 <- t.test(x$jvote[x$tr == 4 & x$ideo7pt == 7], x$jvote[x$tr == 2 & x$ideo7pt == 7])
        
        jf1 <- t.test(x$jfav[x$tr == 4 & x$ideo7pt == 1], x$jfav[x$tr == 2 & x$ideo7pt == 1])
        jf2 <- t.test(x$jfav[x$tr == 4 & x$ideo7pt == 2], x$jfav[x$tr == 2 & x$ideo7pt == 2])
        jf3 <- t.test(x$jfav[x$tr == 4 & x$ideo7pt == 3], x$jfav[x$tr == 2 & x$ideo7pt == 3])
        jf4 <- t.test(x$jfav[x$tr == 4 & x$ideo7pt == 4], x$jfav[x$tr == 2 & x$ideo7pt == 4])
        jf5 <- t.test(x$jfav[x$tr == 4 & x$ideo7pt == 5], x$jfav[x$tr == 2 & x$ideo7pt == 5])
        jf6 <- t.test(x$jfav[x$tr == 4 & x$ideo7pt == 6], x$jfav[x$tr == 2 & x$ideo7pt == 6])
        jf7 <- t.test(x$jfav[x$tr == 4 & x$ideo7pt == 7], x$jfav[x$tr == 2 & x$ideo7pt == 7])
        
        cv1 <- t.test(x$cvote[x$tr == 4 & x$ideo7pt == 1], x$cvote[x$tr == 2 & x$ideo7pt == 1])
        cv2 <- t.test(x$cvote[x$tr == 4 & x$ideo7pt == 2], x$cvote[x$tr == 2 & x$ideo7pt == 2])
        cv3 <- t.test(x$cvote[x$tr == 4 & x$ideo7pt == 3], x$cvote[x$tr == 2 & x$ideo7pt == 3])
        cv4 <- t.test(x$cvote[x$tr == 4 & x$ideo7pt == 4], x$cvote[x$tr == 2 & x$ideo7pt == 4])
        cv5 <- t.test(x$cvote[x$tr == 4 & x$ideo7pt == 5], x$cvote[x$tr == 2 & x$ideo7pt == 5])
        cv6 <- t.test(x$cvote[x$tr == 4 & x$ideo7pt == 6], x$cvote[x$tr == 2 & x$ideo7pt == 6])
        cv7 <- t.test(x$cvote[x$tr == 4 & x$ideo7pt == 7], x$cvote[x$tr == 2 & x$ideo7pt == 7])
        
        cf1 <- t.test(x$cfav[x$tr == 4 & x$ideo7pt == 1], x$cfav[x$tr == 2 & x$ideo7pt == 1])
        cf2 <- t.test(x$cfav[x$tr == 4 & x$ideo7pt == 2], x$cfav[x$tr == 2 & x$ideo7pt == 2])
        cf3 <- t.test(x$cfav[x$tr == 4 & x$ideo7pt == 3], x$cfav[x$tr == 2 & x$ideo7pt == 3])
        cf4 <- t.test(x$cfav[x$tr == 4 & x$ideo7pt == 4], x$cfav[x$tr == 2 & x$ideo7pt == 4])
        cf5 <- t.test(x$cfav[x$tr == 4 & x$ideo7pt == 5], x$cfav[x$tr == 2 & x$ideo7pt == 5])
        cf6 <- t.test(x$cfav[x$tr == 4 & x$ideo7pt == 6], x$cfav[x$tr == 2 & x$ideo7pt == 6])
        cf7 <- t.test(x$cfav[x$tr == 4 & x$ideo7pt == 7], x$cfav[x$tr == 2 & x$ideo7pt == 7])
        
        
        # very lib (row 1)
        TA5Cr1 <- c((jf1$est[1] - jf1$est[2]),
        round(jf1$conf.int, digits = 2),
        (jv1$est[1] - jv1$est[2]),
        round(jv1$conf.int, digits = 2),
        (cf1$est[1] - cf1$est[2]),
        round(cf1$conf.int, digits = 2), 
        (cv1$est[1] - cv1$est[2]),
        round(cv1$conf.int, digits = 2))
        
        
        # lib (row 2)
        TA5Cr2 <- c((jf2$est[1] - jf2$est[2]),
        round(jf2$conf.int, digits = 2),
        (jv2$est[1] - jv2$est[2]),
        round(jv2$conf.int, digits = 2),
        (cf2$est[1] - cf2$est[2]),
        round(cf2$conf.int, digits = 2), 
        (cv2$est[1] - cv2$est[2]),
        round(cv2$conf.int, digits = 2))
        
        
        # somewhat lib (row 3)
        TA5Cr3 <- c((jf3$est[1] - jf3$est[2]),
        round(jf3$conf.int, digits = 2),
        (jv3$est[1] - jv3$est[2]),
        round(jv3$conf.int, digits = 2),
        (cf3$est[1] - cf3$est[2]),
        round(cf3$conf.int, digits = 2), 
        (cv3$est[1] - cv3$est[2]),
        round(cv3$conf.int, digits = 2))
        
        
        # moderate (row 4)
        TA5Cr4 <- c((jf4$est[1] - jf4$est[2]),
        round(jf4$conf.int, digits = 2),
        (jv4$est[1] - jv4$est[2]),
        round(jv4$conf.int, digits = 2),
        (cf4$est[1] - cf4$est[2]),
        round(cf4$conf.int, digits = 2), 
        (cv4$est[1] - cv4$est[2]),
        round(cv4$conf.int, digits = 2))
        
        
        
        # somewhat cons (row 5)
        TA5Cr5 <- c((jf5$est[1] - jf5$est[2]),
        round(jf5$conf.int, digits = 2),
        (jv5$est[1] - jv5$est[2]),
        round(jv5$conf.int, digits = 2),
        (cf5$est[1] - cf5$est[2]),
        round(cf5$conf.int, digits = 2), 
        (cv5$est[1] - cv5$est[2]),
        round(cv5$conf.int, digits = 2))
        
        
        # cons (row 6)
        TA5Cr6 <- c((jf6$est[1] - jf6$est[2]),
        round(jf6$conf.int, digits = 2),
        (jv6$est[1] - jv6$est[2]),
        round(jv6$conf.int, digits = 2),
        (cf6$est[1] - cf6$est[2]),
        round(cf6$conf.int, digits = 2), 
        (cv6$est[1] - cv6$est[2]),
        round(cv6$conf.int, digits = 2))
        
        
        # very cons (row 7)
        TA5Cr7 <- c((jf7$est[1] - jf7$est[2]),
        round(jf7$conf.int, digits = 2),
        (jv7$est[1] - jv7$est[2]),
        round(jv7$conf.int, digits = 2),
        (cf7$est[1] - cf7$est[2]),
        round(cf7$conf.int, digits = 2), 
        (cv7$est[1] - cv7$est[2]),
        round(cv7$conf.int, digits = 2))
        
        
  # here's the table
        TA5C <- rbind(TA5Cr1, TA5Cr2, TA5Cr3, TA5Cr4, TA5Cr5, TA5Cr6, TA5Cr7)
        TA5C <- round(TA5C, digits = 2)
        colnames(TA5C) <- c("J fav est", "J fav CI low", "J fav CI high",
                            "J vote est", "J vote CI low", "J vote CI high", 
                            "C fav est", "C fav CI low", "C fav CI high", 
                            "C vote est", "C vote CI low", "C vote CI high")
  # here's the n. obs column of that table, which I entered by hand, but you can
        # see here: 
        table(x[x$tr == 4 | x$tr == 2,]$ideo7pt)
     
        
           
# APPENDIX 5D: controlling for political knowledge and interest ####

        x20pp$highnewsint <- as.numeric(x20pp$newsint == 1)
        x20pp$hipolknowledge <- as.numeric(x20pp$CC15_310a == 1 & x20pp$CC15_310b == 1)
        x20pp$treat <- as.numeric(x20pp$tr > 2)
        x20pp$ideo_1 <- as.numeric(x20pp$ideo7pt == 1)
        x20pp$ideo_2 <- as.numeric(x20pp$ideo7pt == 2)
        x20pp$ideo_3 <- as.numeric(x20pp$ideo7pt == 3)
        x20pp$ideo_4 <- as.numeric(x20pp$ideo7pt == 4)
        x20pp$ideo_5 <- as.numeric(x20pp$ideo7pt == 5)
        x20pp$ideo_6 <- as.numeric(x20pp$ideo7pt == 6)
        x20pp$ideo_7 <- as.numeric(x20pp$ideo7pt == 7)
        
        lmrobust_jv <- lm(jvote ~ I(treat * ideo_1) + I(treat * ideo_2) + I(treat * ideo_3) + I(treat * ideo_4) + I(treat * ideo_5) + 
                            I(treat * ideo_6) + I(treat * ideo_7) + as.factor(ideo7pt) + highnewsint + hipolknowledge, data = x20pp)
        
        lmrobust_jf <- lm(jfav ~ I(treat * ideo_1) + I(treat * ideo_2) + I(treat * ideo_3) + I(treat * ideo_4) + I(treat * ideo_5) + 
                            I(treat * ideo_6) + I(treat * ideo_7) + as.factor(ideo7pt) + highnewsint + hipolknowledge, data = x20pp)
        
        lmrobust_cv <- lm(cvote ~ I(treat * ideo_1) + I(treat * ideo_2) + I(treat * ideo_3) + I(treat * ideo_4) + I(treat * ideo_5) + 
                            I(treat * ideo_6) + I(treat * ideo_7) + as.factor(ideo7pt) + highnewsint + hipolknowledge, data = x20pp)
        
        lmrobust_cf <- lm(cfav ~ I(treat * ideo_1) + I(treat * ideo_2) + I(treat * ideo_3) + I(treat * ideo_4) + I(treat * ideo_5) + 
                            I(treat * ideo_6) + I(treat * ideo_7) + as.factor(ideo7pt) + highnewsint + hipolknowledge, data = x20pp)
        
        
        
        jv.robustest <- round(summary(lmrobust_jv)$coefficients[,1], digits = 2)
        jv.robustSE.cilow <- round(summary(lmrobust_jv)$coefficients[,1], digits = 2) -
          2*round(summary(lmrobust_jv)$coefficients[,2], digits = 2)
        jv.robustSE.cihigh <- round(summary(lmrobust_jv)$coefficients[,1], digits = 2) +
          2*round(summary(lmrobust_jv)$coefficients[,2], digits = 2)
        jv.robustCI.part1 <- paste("[", round(jv.robustSE.cilow, digits = 2), sep = "")
        jv.robustCI.part2 <- paste(round(jv.robustSE.cihigh, digits = 2), "]", sep = "")
        jv.robustCI <- paste(jv.robustCI.part1, jv.robustCI.part2, sep = ", ")
        
        
        jf.robustest <- round(summary(lmrobust_jf)$coefficients[,1], digits = 2)
        jf.robustSE.cilow <- round(summary(lmrobust_jf)$coefficients[,1], digits = 2) -
          2*round(summary(lmrobust_jf)$coefficients[,2], digits = 2)
        jf.robustSE.cihigh <- round(summary(lmrobust_jf)$coefficients[,1], digits = 2) +
          2*round(summary(lmrobust_jf)$coefficients[,2], digits = 2)
        jf.robustCI.part1 <- paste("[", round(jf.robustSE.cilow, digits = 2), sep = "")
        jf.robustCI.part2 <- paste(round(jf.robustSE.cihigh, digits = 2), "]", sep = "")
        jf.robustCI <- paste(jf.robustCI.part1, jf.robustCI.part2, sep = ", ")
        
        
        cv.robustest <- round(summary(lmrobust_cv)$coefficients[,1], digits = 2)
        cv.robustSE.cilow <- round(summary(lmrobust_cv)$coefficients[,1], digits = 2) -
          2*round(summary(lmrobust_cv)$coefficients[,2], digits = 2)
        cv.robustSE.cihigh <- round(summary(lmrobust_cv)$coefficients[,1], digits = 2) +
          2*round(summary(lmrobust_cv)$coefficients[,2], digits = 2)
        cv.robustCI.part1 <- paste("[", round(cv.robustSE.cilow, digits = 2), sep = "")
        cv.robustCI.part2 <- paste(round(cv.robustSE.cihigh, digits = 2), "]", sep = "")
        cv.robustCI <- paste(cv.robustCI.part1, cv.robustCI.part2, sep = ", ")
        
        
        cf.robustest <- round(summary(lmrobust_cf)$coefficients[,1], digits = 2)
        cf.robustSE.cilow <- round(summary(lmrobust_cf)$coefficients[,1], digits = 2) -
          2*round(summary(lmrobust_cf)$coefficients[,2], digits = 2)
        cf.robustSE.cihigh <- round(summary(lmrobust_cf)$coefficients[,1], digits = 2) +
          2*round(summary(lmrobust_cf)$coefficients[,2], digits = 2)
        cf.robustCI.part1 <- paste("[", round(cf.robustSE.cilow, digits = 2), sep = "")
        cf.robustCI.part2 <- paste(round(cf.robustSE.cihigh, digits = 2), "]", sep = "")
        cf.robustCI <- paste(cf.robustCI.part1, cf.robustCI.part2, sep = ", ")
        
        rownames.robust <- c("intercept", "I(treat * ideo_1)", "I(treat * ideo_2)", "I(treat * ideo_3)", 
                             "I(treat * ideo_4)", "I(treat * ideo_5)", "I(treat * ideo_6)", 
                             "I(treat * ideo_7)", "as.factor(ideo7pt)2",  "as.factor(ideo7pt)3",
                             "as.factor(ideo7pt)4", "as.factor(ideo7pt)5",  "as.factor(ideo7pt)6",
                             "as.factor(ideo7pt)7",
                             "highnewsint", "hipolknowledge")
        
        robust.table <- as.data.frame(cbind(rownames.robust, 
                                            jf.robustest, jf.robustCI, 
                                            jv.robustest, jv.robustCI,
                                            cf.robustest, cf.robustCI, 
                                            cv.robustest, cv.robustCI))
        dv.robust.table <- c("skipme", "Johnson-fav", "Jfav_CI", "Johnson-vote", "Jvote_CI", 
                            "Conley-fav", "Cfav_CI",  "Conley-vote", "Cvote_CI")
        
        robust.table <- rbind(dv.robust.table, robust.table)
        robust.table <- robust.table[,-1]
        
# APPENDIX A10A:  Trustworthiness results ----
    variable <- c("Johnson - Trustworthiness (T-C)", "Conley - Trustworthiness (T-C)", 
                      "Net Trustworthiness (J-C)")
        
    t.test(x$jtrust[x$tr == 1], x$jtrust[x$tr == 3])
    jt <- t.test(x$jtrust[x$tr == 4], x$jtrust[x$tr == 2])
        
    jtest<- jt$estimate[1] - jt$estimate[2]
    jtci.low <- jt$conf.int[1]
    jtci.high <- jt$conf.int[2]
    jtpval <- round(jt$p.value, digits = 2)
    jttrest <- jt$estimate[1]
    jtctest <- jt$estimate[2]
        
        
        
    t.test(x$ctrust[x$tr == 1], x$ctrust[x$tr == 3])
    ct <- t.test(x$ctrust[x$tr == 4], x$ctrust[x$tr == 2])
        
    ctest<- ct$estimate[1] - ct$estimate[2]
    ctci.low <- ct$conf.int[1]
    ctci.high <- ct$conf.int[2]
    ctpval <- round(ct$p.value, digits = 2)
    cttrest <- ct$estimate[1]
    ctctest <- ct$estimate[2]
        
        
        
        
        
    t.test(x$nettrust[x$tr == 1], x$nettrust[x$tr == 3])
    nt <- t.test(x$nettrust[x$tr == 4], x$nettrust[x$tr == 2])
        
    ntest<- nt$estimate[1] - nt$estimate[2]
    ntci.low <- nt$conf.int[1]
    ntci.high <- nt$conf.int[2]
    ntpval <- round(nt$p.value, digits = 2)
    nttrest <- nt$estimate[1]
    ntctest <- nt$estimate[2]
        
        estimate <- c(jtest, ctest, ntest)
        ci.low <- c(jtci.low, ctci.low, ntci.low)
        ci.high <- c(jtci.high, ctci.high, ntci.high)
        pval <- c(jtpval, "0.0000", "0.0000")  # hard code the e-06 and e-11 p.values
        
        mechdata <- as.data.frame(cbind(variable, as.numeric(estimate), as.numeric(ci.low), 
                                        as.numeric(ci.high), as.character(pval)))
        names(mechdata) <- c("variable", "estimate", "ci.low", "ci.high", "pval")
        
        
        arrow_plot_name <- c("Johnson - Trust", "Conley - Trust", "Net Trust (J-C)")
        arrow_means_ct <- c(jtctest, jtctest, ntctest)
        arrow_means_tr <- c(jttrest, cttrest, nttrest)
        mecharrowdata <- as.data.frame(
          cbind(as.character(arrow_means_ct), as.character(arrow_means_tr), 
                as.character(arrow_plot_name)))
        names(mecharrowdata) <- c("arrow_means_ct", "arrow_means_tr", "arrow_plot_name")
        
        mecharrowdata$arrow_means_ct <- as.numeric(as.character(arrow_means_ct))
        mecharrowdata$arrow_means_tr <- as.numeric(as.character(arrow_means_tr))
        
        ## Plot for A10A begins here
        plot(x = NULL, y = NULL, axes = F, 
             # set values of different parameters
             xlim = c(-48, 20),
             ylim = c(1, 3.5),
             xlab = "",
             ylab = "",
             main = "")
        # set margins and letter size
        par(cex=0.8, mai = c(0.5, 0.35, 0.5, 0.35))
        
        # add the 0 vertical line and tick marks
        segments(x0 = 0, y0 = 0, x1 = 0, y1 = 3.5)
        axis(side=1,at=c(-20, -10, 0, 10, 20),tick=TRUE, las=2, cex.axis=0.7)
        # Fill the figure with the information which is inside the 'results' dataframe
        
        # First, add the estimates as points
        points(as.character(mechdata$estimate[3]), 1, type = 'p')
        points(as.character(mechdata$estimate[2]), 2, type = 'p')
        points(as.character(mechdata$estimate[1]), 3, type = 'p')
        
        # Then, add the confidence intervals as lines
        segments(x0=as.numeric(as.character(mechdata$ci.low[3])), y0=1, 
                 x1=as.numeric(as.character(mechdata$ci.high[3])), y1=1)
        segments(x0=as.numeric(as.character(mechdata$ci.low[2])), y0=2, 
                 x1=as.numeric(as.character(mechdata$ci.high[2])), y1=2)
        segments(x0=as.numeric(as.character(mechdata$ci.low[1])), y0=3, 
                 x1=as.numeric(as.character(mechdata$ci.high[1])), y1=3)
        
        # Third, add each variable name
        text(-45, 3, mechdata[1,1], adj = 0, cex=1)
        text(-45, 2, mechdata[2,1], adj = 0, cex=1)
        text(-45, 1, mechdata[3,1], adj = 0, cex=1)
        
        # fourth, add pvalues
        text(18, 3, mechdata$pval[1], cex=1)
        text(18, 2, mechdata$pval[2], cex=1)
        text(18, 1, mechdata$pval[3], cex=1)
        
        # finally, add the subtitles  
        text(18, 3.5, "P Value", font = 2)
        text(-42, 3.5, "Outcome", font = 2)

